Comparison of sequencing data processing pipelines and application to underrepresented African human populations.

BACKGROUND: Population genetic studies of humans make increasing use of high-throughput sequencing in order to capture diversity in an unbiased way. There is an abundance of sequencing technologies, bioinformatic tools and the available genomes are increasing in number. Studies have evaluated and compared some of these technologies and tools, such as the Genome Analysis Toolkit (GATK) and its "Best Practices" bioinformatic pipelines. However, studies often focus on a few genomes of Eurasian origin in order to detect technical issues. We instead surveyed the use of the GATK tools and established a pipeline for processing high coverage full genomes from a diverse set of populations, including Sub-Saharan African groups, in order to reveal challenges from human diversity and stratification. RESULTS: We surveyed 29 studies using high-throughput sequencing data, and compared their strategies for data pre-processing and variant calling. We found that processing of data is very variable across studies and that the GATK "Best Practices" are seldom followed strictly. We then compared three versions of a GATK pipeline, differing in the inclusion of an indel realignment step and with a modification of the base quality score recalibration step. We applied the pipelines on a diverse set of 28 individuals. We compared the pipelines in terms of count of called variants and overlap of the callsets. We found that the pipelines resulted in similar callsets, in particular after callset filtering. We also ran one of the pipelines on a larger dataset of 179 individuals. We noted that including more individuals at the joint genotyping step resulted in different counts of variants. At the individual level, we observed that the average genome coverage was correlated to the number of variants called. CONCLUSIONS: We conclude that applying the GATK "Best Practices" pipeline, including their recommended reference datasets, to underrepresented populations does not lead to a decrease in the number of called variants compared to alternative pipelines. We recommend to aim for coverage of > 30X if identifying most variants is important, and to work with large sample sizes at the variant calling stage, also for underrepresented individuals and populations.

Copy number determination of the gene for the human pancreatic polypeptide receptor NPY4R using read depth analysis and droplet digital PCR.

BACKGROUND: Copy number variation (CNV) plays an important role in human genetic diversity and has been associated with multiple complex disorders. Here we investigate a CNV on chromosome 10q11.22 that spans NPY4R, the gene for the appetite-regulating pancreatic polypeptide receptor Y4. This genomic region has been challenging to map due to multiple repeated elements and its precise organization has not yet been resolved. Previous studies using microarrays were interpreted to show that the most common copy number was 2 per genome. RESULTS: We have investigated 18 individuals from the 1000 Genomes project using the well-established method of read depth analysis and the new droplet digital PCR (ddPCR) method. We find that the most common copy number for NPY4R is 4. The estimated number of copies ranged from three to seven based on read depth analyses with Control-FREEC and CNVnator, and from four to seven based on ddPCR. We suggest that the difference between our results and those published previously can be explained by methodological differences such as reference gene choice, data normalization and method reliability. Three high-quality archaic human genomes (two Neanderthal and one Denisova) display four copies of the NPY4R gene indicating that a duplication occurred prior to the human-Neanderthal/Denisova split. CONCLUSIONS: We conclude that ddPCR is a sensitive and reliable method for CNV determination, that it can be used for read depth calibration in CNV studies based on already available whole-genome sequencing data, and that further investigation of NPY4R copy number variation and its consequences are necessary due to the role of Y4 receptor in food intake regulation.

Abolished InsP3R2 function inhibits sweat secretion in both humans and mice.

There are 3 major sweat-producing glands present in skin; eccrine, apocrine, and apoeccrine glands. Due to the high rate of secretion, eccrine sweating is a vital regulator of body temperature in response to thermal stress in humans; therefore, an inability to sweat (anhidrosis) results in heat intolerance that may cause impaired consciousness and death. Here, we have reported 5 members of a consanguineous family with generalized, isolated anhidrosis, but morphologically normal eccrine sweat glands. Whole-genome analysis identified the presence of a homozygous missense mutation in ITPR2, which encodes the type 2 inositol 1,4,5-trisphosphate receptor (InsP3R2), that was present in all affected family members. We determined that the mutation is localized within the pore forming region of InsP3R2 and abrogates Ca2+ release from the endoplasmic reticulum, which suggests that intracellular Ca2+ release by InsP3R2 in clear cells of the sweat glands is important for eccrine sweat production. Itpr2-/- mice exhibited a marked reduction in sweat secretion, and evaluation of sweat glands from Itpr2-/- animals revealed a decrease in Ca2+ response compared with controls. Together, our data indicate that loss of InsP3R2-mediated Ca2+ release causes isolated anhidrosis in humans and suggest that specific InsP3R inhibitors have the potential to reduce sweat production in hyperhidrosis.

Welander distal myopathy caused by an ancient founder mutation in TIA1 associated with perturbed splicing.

Welander distal myopathy (WDM) is an adult onset autosomal dominant disorder characterized by distal limb weakness, which progresses slowly from the fifth decade. All WDM patients are of Swedish or Finnish descent and share a rare chromosome 2p13 haplotype. We restricted the WDM-associated haplotype followed by whole exome sequencing. Within the conserved haplotype, we identified a single heterozygous mutation c.1150G>A (p.E384K) in T-cell intracellular antigen-1 (TIA1) in all WDM patients investigated (n = 43). The TIA1 protein regulates splicing, and translation through direct interaction with mRNA and the p.E384K mutation is located in the C-terminal Q-rich domain that interacts with the U1-C splicing factor. TIA1 has been shown to prevent skipping of SMN2 exon 7, and we show that WDM patients have increased levels of spliced SMN2 in skeletal muscle cells when compared with controls. Immunostaining of WDM muscle biopsies showed accumulation of TIA1 and stress granulae proteins adjacent to intracellular inclusions, a typical finding in WDM. The combined findings strongly suggest that the TIA1 mutation causes perturbed RNA splicing and cellular stress resulting in WDM. The selection against the mutation is likely to be negligible and the age of the TIA1 founder mutation was calculated to approximately 1,050 years, which coincides with the epoch of early seafaring across the Baltic Sea.

Genome-wide sequencing for the identification of rearrangements associated with Tourette syndrome and obsessive-compulsive disorder.

BACKGROUND: Tourette Syndrome (TS) is a neuropsychiatric disorder in children characterized by motor and verbal tics. Although several genes have been suggested in the etiology of TS, the genetic mechanisms remain poorly understood. METHODS: Using cytogenetics and FISH analysis, we identified an apparently balanced t(6,22)(q16.2;p13) in a male patient with TS and obsessive-compulsive disorder (OCD). In order to map the breakpoints and to identify additional submicroscopic rearrangements, we performed whole genome mate-pair sequencing and CGH-array analysis on DNA from the proband. RESULTS: Sequence and CGH array analysis revealed a 400 kb deletion located 1.3 Mb telomeric of the chromosome 6q breakpoint, which has not been reported in controls. The deletion affects three genes (GPR63, NDUFA4 and KLHL32) and overlaps a region previously found deleted in a girl with autistic features and speech delay. The proband's mother, also a carrier of the translocation, was diagnosed with OCD and shares the deletion. We also describe a further potentially related rearrangement which, while unmapped in Homo sapiens, was consistent with the chimpanzee genome. CONCLUSIONS: We conclude that genome-wide sequencing at relatively low resolution can be used for the identification of submicroscopic rearrangements. We also show that large rearrangements may escape detection using standard analysis of whole genome sequencing data. Our findings further provide a candidate region for TS and OCD on chromosome 6q16.

Genetic adaptation of fatty-acid metabolism: a human-specific haplotype increasing the biosynthesis of long-chain omega-3 and omega-6 fatty acids.

Omega-3 and omega-6 long-chain polyunsaturated fatty acids (LC-PUFAs) are essential for the development and function of the human brain. They can be obtained directly from food, e.g., fish, or synthesized from precursor molecules found in vegetable oils. To determine the importance of genetic variability to fatty-acid biosynthesis, we studied FADS1 and FADS2, which encode rate-limiting enzymes for fatty-acid conversion. We performed genome-wide genotyping (n = 5,652 individuals) and targeted resequencing (n = 960 individuals) of the FADS region in five European population cohorts. We also analyzed available genomic data from human populations, archaic hominins, and more distant primates. Our results show that present-day humans have two common FADS haplotypes-defined by 28 closely linked SNPs across 38.9 kb-that differ dramatically in their ability to generate LC-PUFAs. No independent effects on FADS activity were seen for rare SNPs detected by targeted resequencing. The more efficient, evolutionarily derived haplotype appeared after the lineage split leading to modern humans and Neanderthals and shows evidence of positive selection. This human-specific haplotype increases the efficiency of synthesizing essential long-chain fatty acids from precursors and thereby might have provided an advantage in environments with limited access to dietary LC-PUFAs. In the modern world, this haplotype has been associated with lifestyle-related diseases, such as coronary artery disease.

Characterizing and interpreting genetic variation from personal genome sequencing.

Since the completion of the human genome project, there has been enormous progress in the development of novel technologies for DNA sequencing. The advent of next-generation sequencing technologies now makes it possible to sequence an entire human genome in one or a few experiments. As a consequence, several individual human genomes have now been fully sequenced, using different experimental strategies. Although the protocols differ between the various sequencing technologies, the challenges of analyzing the data, calling variation, and interpreting the results are similar for all platforms. Here, we give an overview of the human genome sequencing projects completed to date. The strategies for aligning sequence reads and extracting information about different types of genetic variation from the sequence data are discussed. Identification of structural variation, such as copy number variation and insertion-deletion variants, can be complex, and there are a plethora of algorithms and analysis tools available. We also give an overview of the challenge of interpreting the whole-genome sequence data both from a technical and clinical perspective.

Characterization of copy number-stable regions in the human genome.

In the past few years the number of copy number variants (CNVs) identified in the human genome has increased significantly, but our understanding of the functional impact of CNVs is still limited. Clinically significant variations cannot easily be distinguished from benign, complicating interpretation of patient data. Multiple studies have focused on analysis of regions that vary in copy number in specific disorders. Here we use the opposite strategy and focus our analysis on regions that never seem to vary in the general population, hypothesizing that these are copy number stable because variations within them are deleterious. Our results show that copy number stable regions are characterized by correlation with a number of genomic features, allowing us to define a list of genomic regions that are dosage sensitive in humans. We find that these dosage-sensitive regions show significant overlap with de novo CNVs identified in patients with intellectual disability or autism. There is also a significant association between copy number stable regions and rare inherited variants in autism patients, but not in controls. Based on this predictive power, we propose that copy number stable regions can be used to complement maps of known CNVs to facilitate interpretation of patient data.

Protein contents in biological membranes can explain abnormal solvation of charged and polar residues.

Transmembrane helices are generally believed to insert into membranes based on their hydrophobicity. Nevertheless, there are important exceptions where polar residues have great functional importance, for instance the S4 helix of voltage-gated ion channels. It has been shown experimentally that insertion can be accomplished by hydrophobic counterbalance, predicting an arginine insertion cost of only 2.5 kcal/mol, compared with 14.9 kcal/mol in cyclohexane. Previous simulations of pure bilayers have produced values close to the pure hydrocarbon, which has lead to spirited discussion about the experimental conditions. Here, we have performed computer simulations of models better mimicking biological membranes by explicitly including protein helices at mass fractions from 15% to 55%, as well as an actual translocon. This has a striking effect on the solvation free energy of arginine. With some polar residues present, the solvation cost comes close to experimental observation at approximately 30% mass fraction, and negligible at 40%. In the presence of a translocon in the membrane, the cost of inserting arginine next to the lateral gate can be as low as 3-5 kcal/mol. The effect is mainly due to the extra helices making it easier to retain hydration water. These results offer a possible explanation for the discrepancy between the in vivo hydrophobicity scale and computer simulations and highlight the importance of the high protein contents in membranes. Although many membrane proteins are stable in pure bilayers, such simplified models might not be sufficiently accurate for insertion of polar or charged residues in biological membranes.

The role of lipid composition for insertion and stabilization of amino acids in membranes.

While most membrane protein helices are clearly hydrophobic, recent experiments have indicated that it is possible to insert marginally hydrophobic helices into bilayers and have suggested apparent in vivo free energies of insertion for charged residues that are low, e.g., a few kcals for arginine. In contrast, a number of biophysical simulation studies have predicted that the bilayer interior is close to a pure hydrophobic environment with large penalties for hydrophilic amino acids--and yet the experimental scales do significantly better at predicting actual membrane proteins from sequence. Here, we have systematically studied the dependence of the free energy profiles on lipid properties, including tail length, saturation, headgroup hydrogen bond strength, and charge, both to see to whether the in vivo insertion can be explained in whole or part from lipid composition of the endoplasmic reticulum (ER) membranes, and if the solvation properties can help interpret how protein function depends on the lipids. We find that lipid charge is important to stabilize charged amino acids inside the bilayer (with implications, e.g., for ion channels), that thicker bilayers have higher solvation costs for hydrophilic side chains, and that headgroup hydrogen bond strength determines how adaptive the lipids are as a hydrophobic/hydrophilic solvent. None of the different free energy profiles are even close to the low apparent in vivo insertion cost, which suggests that regardless of the specific ER membrane composition the current experimental results cannot be explained by normal lipid-type variation.

Titratable amino acid solvation in lipid membranes as a function of protonation state.

Knowledge about the insertion and stabilization of membrane proteins is a key step toward understanding their function and enabling membrane protein design. Transmembrane helices are normally quite hydrophobic so as to efficiently insert into membranes, but there are many exceptions with polar or titratable residues. An obvious example is the S4 helices of voltage-gated ion channels with up to 4 arginines, leading to vivid discussion about whether such helices can insert spontaneously, and if so, what their conformation, protonation state, and cost of insertion really are. To address this question, we have determined geometric and energetic solvation properties for different protonation states of the titrateable amino acids, including hydration, side chain orientation, free energy profiles, and effects on the membrane thickness. As expected, charged states are significantly more expensive to insert (8-16 kcal/mol) than neutral variants (1-3 kcal/mol). Although both sets of values exhibit quite high relative correlation with experimental in vivo hydrophobicity scales, the magnitudes of the in vivo hydrophobicity scales are much lower and strikingly appears as a compressed version of the calculated values. This agrees well with computational studies on longer lipids but results in an obvious paradox: the differences between in vivo insertion and simulations cannot be explained by methodological differences in force fields, possible limited hydrophobic thickness of the endoplasmic reticulum (ER) membrane, or parameters; even anionic lipid head groups (PG) only have limited effect on charged side chains, and virtually none for hydrophobic ones. This leads us to propose a model for in vivo insertion that could reconcile these differences and explain the correlation: if there are considerable hydrophobic barriers inside the translocon, the experimental reference state for the solvation free energy when comparing insertion/translocation in vivo would be quite close to the bilayer environment rather than water.

Position-resolved free energy of solvation for amino acids in lipid membranes from molecular dynamics simulations.

Studies of insertion and interactions of amino acids in lipid membranes are pivotal to our understanding of membrane protein structure and function. Calculating the insertion cost as a function of transmembrane helix sequence is thus an important step towards improved membrane protein prediction and eventually drug design. Here, we present position-dependent free energies of solvation for all amino acid analogs along the membrane normal. The profiles cover the entire region from bulk water to hydrophobic core, and were produced from all-atom molecular dynamics simulations. Experimental differences corresponding to mutations and costs for entire segments match experimental data well, and in addition the profiles provide the spatial resolution currently not available from experiments. Polar side-chains largely maintain their hydration and assume quite ordered conformations, which indicates the solvation cost is mainly entropic. The cost of solvating charged side-chains is not only significantly lower than for implicit solvation models, but also close to experiments, meaning these could well maintain their protonation states inside the membrane. The single notable exception to the experimental agreement is proline, which is quite expensive to introduce in vivo despite its hydrophobicity--a difference possibly explained by kinks making it harder to insert helices in the translocon.

Amino-acid solvation structure in transmembrane helices from molecular dynamics simulations.

Understanding the solvation of amino acids in biomembranes is an important step to better explain membrane protein folding. Several experimental studies have shown that polar residues are both common and important in transmembrane segments, which means they have to be solvated in the hydrophobic membrane, at least until helices have aggregated to form integral proteins. In this work, we have used computer simulations to unravel these interactions on the atomic level, and classify intramembrane solvation properties of amino acids. Simulations have been performed for systematic mutations in poly-Leu helices, including not only each amino acid type, but also every z-position in a model helix. Interestingly, many polar or charged residues do not desolvate completely, but rather retain hydration by snorkeling or pulling in water/headgroups--even to the extent where many of them exist in a microscopic polar environment, with hydration levels corresponding well to experimental hydrophobicity scales. This suggests that even for polar/charged residues a large part of solvation cost is due to entropy, not enthalpy loss. Both hydration level and hydrogen bonding exhibit clear position-dependence. Basic side chains cause much less membrane distortion than acidic, since they are able to form hydrogen bonds with carbonyl groups instead of water or headgroups. This preference is supported by sequence statistics, where basic residues have increased relative occurrence at carbonyl z-coordinates. Snorkeling effects and N-/C-terminal orientation bias are directly observed, which significantly reduces the effective thickness of the hydrophobic core. Aromatic side chains intercalate efficiently with lipid chains (improving Trp/Tyr anchoring to the interface) and Ser/Thr residues are stabilized by hydroxyl groups sharing hydrogen bonds to backbone oxygens.